Solutions Day 4

welcome to class!

Data Introduction

Before we begin, let’s bring in the data we’ve been working with.

df <- read_csv("https://www.fema.gov/api/open/v2/DisasterDeclarationsSummaries.csv")

df_new <- df |> 
  mutate(GEOID=str_c(fipsStateCode, fipsCountyCode))

county_pop <- read_csv("data-raw/county_population.csv")

joined_new <- left_join(df_new, county_pop, by="GEOID") |> 
  mutate(year=year(incidentBeginDate))

Let’s take a look at what were working with. Check out joined_new with the usual function we use:

glimpse(joined_new)

Rows: 64,955
Columns: 31
$ femaDeclarationString    <chr> "FM-5389-AZ", "FM-5389-AZ", "FM-5464-RI", "FM…
$ disasterNumber           <dbl> 5389, 5389, 5464, 5463, 5462, 4731, 5460, 545…
$ state                    <chr> "AZ", "AZ", "RI", "KS", "NE", "CO", "OK", "OK…
$ declarationType          <chr> "FM", "FM", "FM", "FM", "FM", "DR", "FM", "FM…
$ declarationDate          <dttm> 2021-06-06, 2021-06-06, 2023-04-14, 2023-04-…
$ fyDeclared               <dbl> 2021, 2021, 2023, 2023, 2023, 2023, 2023, 202…
$ incidentType             <chr> "Fire", "Fire", "Fire", "Fire", "Fire", "Floo…
$ declarationTitle         <chr> "TELEGRAPH FIRE", "TELEGRAPH FIRE", "QUEENS R…
$ ihProgramDeclared        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ iaProgramDeclared        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ paProgramDeclared        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ hmProgramDeclared        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ incidentBeginDate        <dttm> 2021-06-06, 2021-06-06, 2023-04-14, 2023-04-…
$ incidentEndDate          <dttm> NA, NA, 2023-04-16, 2023-04-16, NA, 2023-06-…
$ disasterCloseoutDate     <dttm> 2023-09-29, 2023-09-29, NA, NA, NA, NA, NA, …
$ tribalRequest            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ fipsStateCode            <chr> "04", "04", "44", "20", "31", "08", "40", "40…
$ fipsCountyCode           <chr> "007", "021", "009", "201", "025", "009", "14…
$ placeCode                <dbl> 99007, 99021, 99009, 99201, 99025, 99009, 991…
$ designatedArea           <chr> "Gila (County)", "Pinal (County)", "Washingto…
$ declarationRequestNumber <dbl> 21041, 21041, 23042, 23038, 23036, 23081, 230…
$ lastIAFilingDate         <dttm> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ lastRefresh              <dttm> 2023-09-29 21:02:14, 2023-09-29 21:02:14, 20…
$ hash                     <chr> "14b4a2c314124cae33e8ec790782a12c2af10a4c", "…
$ id                       <chr> "226d44a3-418a-4102-9a55-385ea6599b2a", "7304…
$ GEOID                    <chr> "04007", "04021", "44009", "20201", "31025", …
$ NAME                     <chr> "Gila County, Arizona", "Pinal County, Arizon…
$ variable                 <chr> "B01003_001", "B01003_001", "B01003_001", "B0…
$ estimate                 <dbl> 53846, 447559, 126139, 5474, 26022, 3570, 519…
$ moe                      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, 164, NA, …
$ year                     <dbl> 2021, 2021, 2023, 2023, 2023, 2023, 2023, 202…

Okay, let’s transform the data like we did before.

Can you count up how many disasters there have been per year? Not total disasters. Individual disasters.

Call the new column “total”.

annual_disasters <- joined_new |> 
  count(incidentType, year, name="total") 
  
    
annual_disasters

You only need to add one new line. 
The function starts with a *c* and don't forget to name the column you're creating.
year should be the second argument in the function.

Alright, we’ve got a lot of data going back decades.

The benefit of working quickly with data in R is how you can quickly visualize it to spot any trends.

Let’s do that.

But before we do, let’s create another data frame specifically for fires.

Filter incidentType for “Fire”, please.

annual_fires <- annual_disasters  |> 
  filter(incidentType=="Fire")

annual_fires

The function starts with a *f* and don't forget ==

Grammar of Graphics

The grammar of graphics lets you approach visualizations structurally, letting you combine and swap out graphical elements into figures that display data meaningfully.

It takes two lines of code.

This is what the code and chart looks like.

Run the code below.

ggplot(data=annual_fires) +
  geom_col(aes(x=year, y=total)) 

Basically, every of a chart can created using these aesthetic components and mapping them:

Mapping data to aesthetics

Aesthetic

• The visual property of a graph

• Position, shape, color, etc.

Data

• A column in a data set

---

Here’s are the core components of the chart above:

Data	Aesthetic	Graphic/Geometry
Year	Position (x-axis)	Column
Total disasters	Position (y-axis)	Point

---

Here’s how the data was mapped in ggplot2 code from the annual_fires data frame:

Data	aes()	geom
year	x	geom_col()
total	y	geom_col()

ggplot() template

Here’s the dataframe called annual_fires as a reminder:

annual_fires |> slice(1:5)

Okay, now that you see where all the pieces come from, here’s how ggplot() works.

At its core you need to tell it what data you’re using, what type of visual geometry you want to use, and what variables you want represented from the data.

Important: We have to use + signs between each line, not |>. This is because ggplot() was created before the tidyverse piping method was established.

---

Grammatical layers

When constructing charts, so far we know about data, aesthetics, and geometries.

Think of these components as layers.

Add them to foundational ggplot() with +

These are all the arguments we can enhance the data viz with.

Change the colors of the viz based on a column. Or the size of the shape.

Or the opacity or the gradient.

Possible aesthetics

We can also swap out the different geometry types.

If you don’t want a bar, you can use a line. Or a point.

You can even use shapefiles.

Possible geoms

THERE ARE SO MANY GEOMS for different visualizations. Here are the official ones.

Try the code from above again but this time use geom_point() and then try it with geom_line()

ggplot(data=annual_fires) +
  geom_point(aes(x=year, y=total)) 

You can really start to see the power of cycling quickly through different chart styles to see which one is most effective at telling the story you want to tell.

So after you have the very basic elements needed to create a chart, you can build and style it with more layers.

Because the defaults are rarely what you want and effective dataviz comes from small decisions you make along the way.

Additional layers

There are many of other grammatical layers we can use to describe graphs.

We sequentially add layers onto the foundational ggplot() plot to create complex figures.

Scales change the properties of the variable mapping.

Here are a few examples:

Example layer	What it does
scale_x_continuous()	Make the x-axis continuous
scale_x_continuous(breaks = 1:5)	Manually specify axis ticks
scale_x_date()	Considers x-axis dates
scale_color_gradient()	Use a gradient
scale_fill_viridis_d()	Fill with discrete viridis colors

---

Check out the x-axis.

Exercise 2

Now add scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) to the bottom of the code.

ggplot(data=annual_fires) +
  geom_col(aes(x=year, y=total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022)

Do you see the difference at the bottom of the chart compared to the one above it?

It limited the scope of the x-axis so it didn’t go back to the ’50s anymore.

And it specifically labeled the years 2010 through 2022.

Facets

The next possible layer allows for small multiples. It’s really neat.

Facets show subplots for different subsets of data.

Example layer	What it does
facet_wrap(vars(incidentType))	Plot for each disaster type
facet_wrap(vars(incidentType, year))	Plot for each disaster type/year
facet_wrap(…, ncol = 1)	Put all facets in one column
facet_wrap(…, nrow = 1)	Put all facets in one row

The table above shows all the different ways you can use facets– you can break it out by one extra variable or even two.

We’ll use the annual disasters this time so we have more than just the fires.

But we’ll filter it to hurricanes and fires and floods.

And we can combine it with pipes before we use ggplot() it.

Add the facet_wrap() line on the variable incidentType (like the first example in the table above).

annual_disasters |> 
  filter(incidentType %in% c("Hurricane", "Fire", "Flood")) |> 
ggplot() +
  geom_col(mapping=aes(x= year, y= total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) +
  facet_wrap(vars(incidentType))

Alright, looks like the x-axis labels are getting a little crowded.

We can’t even read it!

Try again!

Now, try it with ncol=1 as an additional argument in facet_wrap()

annual_disasters |> 
  filter(incidentType %in% c("Hurricane", "Fire", "Flood")) |> 
ggplot() +
  geom_col(mapping=aes(x= year, y= total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) +
  facet_wrap(vars(incidentType), ncol=1)

Labels

Now we can add more customization to the chart.

To make it really shine!

Example layer	What it does
labs(title = “Neat title”)	Title
labs(caption = “Something”)	Caption
labs(y = “Something”)	y-axis
labs(color = “Type”)	Title of size legend

• Title should be “Disaster declarations since 2010”
• Label for the x-axis should be blank (aka ““) because the years are obvious
• Label for the y-axis should be “Total”
• Caption should be “Data: FEMA”

Add those labels below:

annual_disasters |> 
  filter(incidentType %in% c("Hurricane", "Fire", "Flood")) |> 
ggplot() +
  geom_col(mapping=aes(x= year, y= total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) +
  facet_wrap(vars(incidentType), ncol=1) +
  labs(
    title = "Disaster declarations since 2010",
    x = "",
    y = "Total",
    caption= "Data: FEMA"
  )

You only need to call labs() once.
Within parentheses, just separate the arguments with commas. You don't use the plus signs.

Themes

Change the appearance of anything in the plot.

While you can customize every font, color, gradient, etc, you can set these styles up ahead of time or use the ones others have created.

There are many built-in themes.

Example layer	What it does
theme_grey()	Default grey background
theme_bw()	Black and white
theme_dark()	Dark
theme_minimal()	Minimal

Try out the different themes listed above in the code below.

annual_disasters |> 
  filter(incidentType %in% c("Hurricane", "Fire", "Flood")) |> 
ggplot() +
  geom_col(mapping=aes(x= year, y= total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) +
  facet_wrap(vars(incidentType), ncol=1) +
  labs(
    title = "Disaster declarations since 2010",
    x = "",
    y = "Total",
    caption= "Data: FEMA"
  ) +
  theme_minimal()

More themes

There are a collections of pre-built themes online, like the ggthemes package.

Organizations often make their own custom themes, like the BBC.

Theme adjustments

Make theme adjustments with theme()

There are a billion options here!

Add this chunk of code in the exercise below it:

theme_bw() + 
theme(plot.title = element_text(face = "bold"),
      panel.grid = element_blank(),
      axis.title.y = element_text(face = "italic"))

Exercise 5

annual_disasters |> 
  filter(incidentType %in% c("Hurricane", "Fire", "Flood")) |> 
ggplot() +
  geom_col(mapping=aes(x= year, y= total)) +
  scale_x_continuous(limits=c(2010, 2022), breaks=2010:2022) +
  facet_wrap(vars(incidentType), ncol=1) +
  labs(
    title = "Disaster declarations since 2010",
    x = "",
    y = "Total",
    caption= "Data: FEMA"
  ) +
theme_bw() + 
theme(plot.title = element_text(face = "bold"),
      panel.grid = element_blank(),
      axis.title.y = element_text(face = "italic"))

Warning: Removed 144 rows containing missing values (`position_stack()`).

Warning: Removed 6 rows containing missing values (`geom_col()`).

These were just a few examples of layers.

See the ggplot2 documentation for complete examples of everything you can do

Done!

Congrats on completing the walkthroughs for Class 4!

But there’s more

Putting it all together

Usually we’d next go over all the different geom_ visualizations you can create using a single data set.

• A ggplot2 tutorial for beautiful plotting in R

But we’re going to use more real-life data that I think will be relevant to your journalism.

The data set is raw deaths data from San Diego. It’s a combination of 1997-2019 data from San Diego’s data portal and 2020 data from a public information request on MuckRock.

Downloads this data san_diego.csv and place it in your project folder.

Make sure you’ve got the proper libraries loaded.

Ready for the code?

library(tidyverse)
library(janitor)
library(lubridate)

sd <- read_csv("data-raw/san_diego.csv")

Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)

Rows: 66218 Columns: 26
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (21): Death Date, Security Status, Gender, Race, Ethnic Group (Standardi...
dbl  (5): Year, Age in Years, Event Zip, Death Zip, Res Zip

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

sd <- clean_names(sd)

sd_adjusted <- sd |> 
  mutate(death_date=mdy(death_date)) |> 
  mutate(month=month(death_date, label=TRUE, abbr=TRUE)) 

sd_adjusted

Alright, I’ve cleaned it up for you.

There’s some really great data here. It’s got gender, race, and several levels of manner of death stretching back to 1997 and through possibly October of 2020.

Let’s start summarizing the data so we can start looking for trends.

Exercise 7

Can you count up the number of deaths by manner_of_death by month and year, please?

sd_month <- sd_adjusted |> 
  count(year, month, manner_of_death, name="deaths") |> 
  mutate(date=mdy(paste0(month, " 1, ", year)))

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `date = mdy(paste0(month, " 1, ", year))`.
Caused by warning:
!  35 failed to parse.

sd_month

Hint: count / deaths

Exercise 8

Now, fill in the blank below to create a line chart for each type of death

sd_month |> 
  ggplot(aes(x=date, y=deaths, color=manner_of_death, group=manner_of_death)) +
  geom_line() +
  labs(title="Monthly deaths in San Diego",
       caption="Source; San Diego Medical Examiner")

Warning: Removed 35 rows containing missing values (`geom_line()`).

Alright, a little messy. We can see some overall growing trend in one category, but that’s it.

Perhaps looking at the data by month is too granular. Let’s step back and aggregate by year.

By year

I went ahead and created a new sd_year dataframe counting up the deaths by year (while excluding October, November, and December) so we can compare prior years to this year.

sd_year <- sd_adjusted |> 
  # if we're going to compare this year to previous years, we need to exclude data we don't have yet
  filter(!month %in% c("Oct", "Nov", "Dec")) |> 
  count(year, manner_of_death, name="deaths") 

datatable(sd_year)

Exercise 9

Okay, your turn to make a chart.

Make me a faceted chart that breaks out all the individual manner_of_death into its own chart, small-multiple style.

sd_year |> ggplot(aes(x=year, y=deaths)) +
  geom_col() +
  facet_wrap(vars(manner_of_death), ncol=4)

Alright, now we’re getting somewhere.

Looks like accidents have been trending up year over year.

If we focus on 2020, it looks like Natural causes have increased. But it also increased in 2019.

Suicides actually look down this year compared to prior years.

Hm…

What else can we do?

We can try to measure Excess Deaths

Average each month by every year prior to 2020 and compare it to 2020’s trend line.

I’ll give you the code again.

We’re going to use a function called case_when to create a new column called year_type. If the year is 2020, then it will be “2020” otherwise it will be “1997-2020”. And then we find the average number of deaths for each month for those two groups.

sd_group <- sd_adjusted |> 
  filter(!month %in% c("Oct", "Nov", "Dec")) |> 
  count(year, month, manner_of_death, name="deaths") |> 
  mutate(year_type=case_when(
    year==2020 ~ "2020",
    TRUE ~ "1997-2019"
  )) |> 
  group_by(month, manner_of_death, year_type) |> 
  summarize(avg_deaths=mean(deaths, na.rm=T)) |> 
  filter(!is.na(month))

datatable(sd_group)

Looking very smooth.

Let’s chart it.

Exercise 10

Can you please create a faceted line chart of the data above? But with year_type as two different lines?

Fill in the three blanks to generate the chart.

ggplot(sd_group, aes(x=month, y=avg_deaths, color=year_type, group=year_type)) +
  geom_line() +
  facet_wrap(vars(manner_of_death), scales="free_y", ncol=2)

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

Great.

Now there’s some interesting trend in July and August, right?

And it shows that maybe the last month of data is so low. It’s possible that the data for the month is incomplete and should be excluded from analysis.

How about we come up with a number for the nutgraf of a story?

How many more accidental cause deaths are there in 2020 compared to the historical average?

Exclude September since it seems so off

Here’s a look at the dataframe you can estimate this with.

datatable(sd_group)

What code do you put below to prove this sentence:

“There were X% more accidental deaths in in Jan 2020 than the historical average of January in San Diego”

hint: you’ll need filter(), pivot_wider(), group_by(), summarize(), mutate(), and math

sd_group  |> 
  filter(month=="Jan" & manner_of_death=="Accident") |>
  pivot_wider(names_from="year_type", values_from="avg_deaths") %>% 
  mutate(change=(`2020`-`1997-2019`)/`1997-2019`*100)

Good job!

Further exploratory visual analysis

Alright, comparing 2020 to the average from previous years seems to be a good decision.

Some interesting trends that we could actually write about are surfacing.

Let’s compare that overall instead of by month.

Here’s the code. We’re also excluding September based on what we discovered in the chart above.

sd_group_compare <- sd_adjusted |> 
  filter(!month %in% c("Sep", "Oct", "Nov", "Dec")) |> 
  count(year, manner_of_death, name="deaths") |> 
  mutate(year_type=case_when(
    year==2020 ~ "2020",
    TRUE ~ "1997-2019"
  )) |> 
  group_by(manner_of_death, year_type) |> 
  summarize(avg_deaths=round(mean(deaths, na.rm=T)))

`summarise()` has grouped output by 'manner_of_death'. You can override using
the `.groups` argument.

datatable(sd_group_compare)

Exercise 11

Run the code below with manner_of_death as x and avg_deaths as y.

Then swap them.

Which do you prefer and why?

ggplot(sd_group_compare, aes(x=manner_of_death, y=avg_deaths, fill=year_type)) +
  geom_bar(position="dodge", stat="identity") 

Alright, before we go, I want to clean things up some.

I want to get rid of the manners of death that have barely any and I want to reorder the labels so that it’s in alphabetical order.

Exercise 12

Take a look at the code below. Absorb it.

Then generate the code and see what pops up.

sd_group_compare |> 
  filter(!manner_of_death %in% c("Other", "Family Paid Autopsy")) |> 
  filter(!is.na(manner_of_death)) |> 
  ggplot(aes(x=avg_deaths, y=forcats::fct_rev(manner_of_death),  fill=year_type)) +
  geom_bar(position="dodge", stat="identity") +
  labs(title="Manner of death in San Diego",
       subtitle="January and August deaths in 2020 compared to average deaths between 1997 and 2019",
       caption="Source: San Diego Medical Examiner",
       y="Manner of death",
       x="Average deaths",
       fill="Year") +
  theme_minimal()

Story

So, what do you think the story is?

In San Diego, accidents are way up, suicides are slightly up, and meanwhile homicides are down.

What can we do next?

Well, dig into the accidents, perhaps, and see if there’s any explanation for the huge increase.

Alright, congratulations on going on this exploratory data visualization journey.

Some of the answers won’t appear right away unless you poke around and look at the data in as many ways as possible.